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ABSTRACT 

In  this  paper  we  present  a  class  of  multiresolution  algorithms  for  fast  application  of 
structured  dense  matrices  to  arbitrary  vectors,  which  includes  the  fast  wavelet  transform 
of  Beylkin,  Coifinan  and  Rokhlin  and  the  multilevel  matrix  multiplication  of  Brandt  and 
Lubrecht.  In  designing  these  algorithms  we  first  apply  data  compression  techniques  to  the 
matrix  and  then  show  how  to  compute  the  desired  matrix-vector  multiplication  from  the 
compressed  form  of  the  matrix.  In  describing  this  class  we  pay  special  attention  to  an 
algorithm  which  is  based  on  discretization  by  cell-averages  as  it  seems  to  be  suitable  for 
discretization  of  integral  transforms  with  integrably  singular  kernels. 
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1.  Introduction. 


In  this  paper  we  present  a  class  of  multiresolution  algorithms  for  rapid  application 
of  dense  matrices  to  vectors.  A  direct  application  of  an  arbitrary  NxN  dense  matrix 
to  a  vector  requires  operations.  However,  when  the  matrix-vector  multiplication 
stems  from  a  discretization  of  an  integral  transform 

(1.1)  u(t)  =  J  J  K(x,y)v(y)dy, 

where  the  kernel  K{x,y)  is  smooth  except  possibly  along  curves,  this  product  can 
be  performed  to  any  prescribed  accuracy  with  only  0{N)  operations. 

In  [1]  Beylkin,  Coifman  and  Rokhlin  (BCR)  present  a  wavelet  based  algorithm 
(referred  to  as  the  “nonstandard  form”),  in  which  the  matrix- vector  multiplication 
is  performed  by  successive  contributions  from  different  scales.  It  starts  with  an 
initial  blurred  (low  resolution)  output  vector  for  u  in  (1.1),  which  is  then  upgraded 
successively  to  higher  and  higher  resolution,  in  much  the  same  way  as  the  pyramid 
scheme  in  image  comrpession. 

In  [2]  Brandt  and  Lubrecht  (BL)  describe  a  multilevel  matrix-vector  multiplica¬ 
tion  which  is  viewed  as  performing  part  of  the  integration  in  (1.1)  on  coeirser  grids. 
This  is  possible  wherever  the  local  smoothness  of  the  kernel  K{x,  y)  enables  the  re¬ 
placement  of  its  fine  grid  values  by  sufficiently  accurate  interpolation  from  coarser 
grids. 

In  [7]  we  have  presented  a  class  of  multiresolution  edgorithms  for  data  compres¬ 
sion.  In  the  present  paper  we  apply  these  data  compression  edgorithms  to  matrices 
as  a  tensor  product  of  one-dimensional  oeprators  to  obtain  a  multiresolution  rep¬ 
resentation  of  the  matrix.  Using  this  representation  we  derive  a  class  of  0{N) 
matrix-vector  multiplication  algorithms,  which  includes  the  BCR  algorithm  [1]  and 
the  BL  algorithm  [2]  as  paurticuleu:  cases.  In  describing  this  class  we  also  pay  specied 
attention  to  the  edgorithm  which  is  based  on  discretization  by  cell- averages,  because 
it  seems  to  be  particularly  suitable  to  kernels  with  integrable  singularity. 

2.  Discretization  and  Reconstruction. 

In  this  section  we  describe  a  class  of  discretizations  of  a  function  and  the  approxi¬ 
mate  inverse  of  these  discretizations,  namely  the  approximate  recovery  of  a  function 
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from  its  given  discrete  values;  we  refer  to  the  process  of  recovery  as  reconstruction. 

Let  } ,  Xj  =  j-k°  be  a  partition  of  the  real  line  into  uniform  intervals  {/j  } ,  Ij  = 
of  size  ho.  Let  ip(x)  be  a  function  which  is  concentrated  around  z  =  0 
and  satisfies 

(2.1a) 

and  define  its  scaled  translates 

(2.1b)  = 

Given  a  function  f(x)  we  discretize  it  by 

(2.2)  ^  =  (/,v>“)  =  J 

Next  let  us  introduce  an  approximate  recovery  of  the  function  /(z)  from  its  given 
values  /°  =  {fj}  which  we  refer  to  as  reconstruction  and  denote  by  7Z.{x;f^).  We 
say  that  the  reconstruction  is  r-th  order  accurate  if 

(2.3a)  7l{x]p)  =  /(z)  +  O((ho)’’),  (accuracy) 

provided  that  /(z)  is  sufficiently  smooth.  We  assume  that  the  reconstruction  is 
conservative  in  the  sense  that 

(2.3b)  {1Z{-;p),  (fj)  =  ^  (conservation). 

Unlike  the  setup  in  [7],  where  the  goal  is  to  obtain  maximal  data  compression,  in 
the  present  application  to  matrix- vector  multiplication  we  want  minimal  number  of 
operations.  Therefore  we  assume  that  'R.{-\p)  is  a  linear  functional  of  /°  and  that 
(^(z)  satisfies  a  dilation  equation 

(2.4a)  (p{x)  =  2  ^  a/(^(2z  —  £), 

t 
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where  the  coefficients  {a/}  satisfy 


(2.4b) 

(2.4c) 


1 

aioii+2m  =0  for  m  ^  0. 


We  note  that  relation  (2.4b)  is  just  a  consistency  condition.  Given  a  set  of  {a/} , 

1,  it  is  shown  in  [4]  and  [8]  that  <^{x)  is  determined  by  the  dilation  equation  (2.4a) 
up  to  a  multiplicative  constant.  Hence  (p{x)  is  determined  uniquely  by  adding  the 
normalization  (2.1a)  to  (2.4a)-(2.4b).  In  Appendix  A  we  show  that  condition  (2.4c) 
implies  orthogonality  of  some  matrices  and  thus  reduces  the  number  of  operations 
in  our  algorithm.  In  order  for  the  set  of  functions  {<^5(3^)}  to  be  orthogonal  we  have 
to  add  zuiother  consistency  relation  (see  [8]) 


(2.5a) 


in  which  case 


(2.5b) 


/  0  Ov 


where  Sij  is  the  Kronecker-^;  i.e.  6i^i  =  1,  Sij  =  0  for  i  ^  j. 

In  this  paper  we  highlight  the  following  three  cases: 

Case  1.  Point  values. 


(2.6a) 


ip{x)  =  S(x) 


where  ^(a:)  is  Dirac’s  distribution.  As  pointed  out  by  Strang  [8]  it  satisfies  the 
dilation  equation 


(2.6b) 


S(x)  =  2S(2x) 


and  thus 


(2.6c) 


ao  =  1,  =  0  for  £  ^  0. 
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Note  that  the  coefficients  (2.6c)  trivially  satisfy  the  orthogonality  relation  (2.4c). 
However 

E“?  =  i 

and  thus  (2.5)  is  not  valid  in  this  case. 

The  discretization  (2.2)  becomes 

(2.7a) 

i.e.  the  function  f(x)  is  discretized  by  taking  its  value  at  the  grid  points  {x°}.  The 
conservation  property  (2.3b)  becomes 

(2-76)  7l(x;;/°)  =  /7, 

i.e.  the  reconstruction  is  an  interpolation  of  the  values  {/“}  at  the  grid  points  {zj}. 
Case  2.  CelKaveiages. 

r  1  -1  <z  <  0 

(2.8a)  ip{x)  =  A'[_i,o)(x)  = 

1.  0  otherwise, 

satisfies  the  dilation  equation 

(2.8b)  ip(x)  =  <^(2z)  +  (p{2x  +  1) 

and  thus 

(2-8c)  ao  =  a_i  =  ^,  a/ =  0  for  £  ^  — 1,0. 

The  discretization  of  /(z)  in  (2.2)  becomes 

(2.9a)  /7  =  (±  -,)  ^ 
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i.e.  f{x)  is  discretized  by  taking  /°  to  be  its  average  in  the  interval  7®.  The 
conservation  requirement  (2.3b)  becomes 


(2.9b) 


-Rix-,r)dz=rj 


Let  us  denote  by  F{x) 


(2.10a) 


F{x)  =  r 

Jo 


the  primitive  function  of  f(x) 


(2.10b) 


£f(i)  =  f(x) 


and  observe  that 


(2.10c) 


It  is  easy  to  see  that 


(2.11) 


Tl(x-,p)  =  ^I{x-,F^), 


where  7(x;F®)  is  any  interpolation  of  the  values  F®  =  F(xj)  (2.10c),  satisfies 
the  conservation  requirement  (2.9b).  This  reconstruction  procedure  is  r-th  order 
accurate  (2.3a)  if  the  interpolation  technique  in  (2.11)  satisfies 


(2.12) 


^/(i;  F»)  =  £f(x)  +  O((hor)  =  /(x)  +  O((hor) 


for  sufficiently  smooth  f(x). 

Case  S.  Orthogonal  Wavelets. 

Let  <^(x)  be  a  function  which  is  determined  by  the  dilation  equation  (2.4a),  with 
coefficients  that  satisfy  (2.4b)-(2.4c)  and  (2.5a).  Thus  we  assume  orthogonality  of 
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the  set  {(^°}  (2.5b).  In  the  context  of  this  paper  it  is  most  natural  to  describe 
wavelets  by  first  specifying  the  reconstruction  to  be  a  linear  combination  of 
i.e. 


(2.13a)  = 

i 

and  to  leave  the  discretization  (2.2)  to  be  determined  later.  The  conservation  re¬ 
quirement  (2.3b)  becomes 

(2.13b)  (n-,P),  v°i)  =  =  Pi- 

i 

Using  the  orthogonality  (2.5b)  we  get 
(2.13c)  a.  = 

Thus 

(2.14) 

Using  the  theory  of  approximation  by  translates  Strang  [8]  shows  that  in  order  for 
the  reconstruction  (2.14)  to  be  r-th  order  accurate  (2.3a)  we  have  to  impose  the 
following  condition  on  the  coefficients  {a/}, 

(2.15)  ^(-l)'^’"a/ =  0  form  =  0,l,...  ,r-l. 

Daubechies  [4]  showed  that  in  order  to  satisfy  the  conditions  on  {o/}  listed  so 
far,  one  needs  at  least  2r  nonzero  coefficients,  and  that  the  set  of  exactly  2r  nonzero 
coefficients  is  unique.  For  r  =  1  this  solution  is  given  by  (2.8c),  i.e.  v’(z)  is  the 
box  function  (2.8a).  For  r  >  2  the  resulting  ip{x)  is  necessarily  nonsymmetric; 
the  smoothness  of  <^(z)  increjises  with  r,  but  only  by  half  a  derivative  (approxi¬ 
mately)  each  time.  Beylkin,  Coifman  and  Rokhlin  in  [1]  impose  an  additional  set 
of  requirements  on  {o/},  namely  that  there  exists  an  integer  so  that 

(2.16a)  f  ip{x  +  Tr)x'^ dx  =  0  for  m  =  1,2, . . .  ,r  -  1; 
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This  implies 

(2.16b)  =  (/,  ^5)  =  /(i;  +  T.ho) + o((h,n 

which  shows  that  the  integration  in  (2.16b)  can  be  approximated  to  r-th  order 
accuracy  by  a  single  point  quadrature.  They  show  that  there  is  a  solution  to  the 
extended  set  of  conditions  with  3r  nonzero  coefficients  {oct}. 

We  remark  that  for  large  r,  the  discretization  implied  by  (2.16b)  is  close  to  that 
of  point  values. 

3.  Multiresolution  Algorithms  for  Data  Compression 

In  this  section  we  consider  a  situation  where  we  are  given  No  values 

(3.1a)  P  =  Nq=2^\  no  integer, 

which  represent  a  discretization  (2.2)  of  some  function  f{x)  corresponding  to  a 
uniform  partition  of  [0, 1], 

(3.1b)  Xj=j'  ho,  0  <  j  <  No,  ho  =  1/A^o- 

To  simplify  our  presentation  we  2issume  for  the  time  being  that  /(i)  is  periodic 
with  period  1,  so  that  values  outside  [0,1]  are  known  by  periodic  extension. 

We  consider  the  set  of  nested  grids 

(3.2a)  ht  =  l/Nt,  ^It  =  2-‘JV„; 

for  0  <  <  X,  where  k  =  0,  the  original  grid,  is  the  finest  in  the  hierarchy  and 

k  —  L,  L  <  no,  is  the  coarsest.  The  coarser  (k  +  l)-th  grid  is  formed  form  the  k-th 
grid  by  removing  the  grid  points  {ijj-i  thus 

(3.2b)  0<j<Nk+u  Nk^,=Nk/2. 

To  each  of  the  nested  grids  we  associate  a  discretization 

(3.3a)  =  //  =  (/.¥'*), 
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where  (pj  is  properly  scaled 


(3.3b) 


It  follows  from  the  dilation  relation  (2.4a)  that 


(3.3c) 


and  consequently 


(3.3d) 


i)  =  1  s  J  < 


We  rewrite  (3.3d)  in  the  matrix  form 


/*  =  Hn,  X  2N.. 


Given  /°  we  use  (3.4)  to  successively  compute  P,...  ,/^.  Observe  that  these 
values  are  not  computed  from  the  definition  (3.3a)  but  from  the  dilation  relation 
(3.3d);  thus  no  explicit  knowledge  of  f{x)  is  required. 

Given  /*  we  can  use  the  reconstruction  ‘7J(x;  /* )  in  order  to  get  an  approximation 
to  the  discrete  values  p~^  of  the  finer  level  by 


(3.5a) 


ft'  =  ’  <  J  <  2iVn  =  JVi-,. 


As  we  have  mentioned  earlier,  in  this  paper  we  take  ??.(•;  /)  to  be  a  linear  functional 
of  /.  Hence  (3.5a)  can  be  epxressed  in  the  matrix  form 


(3.5b) 


/*-*  =Rp 


where  A  is  an  2Nk  x  Nk  matrix.  Because  of  (3.3c)  and  the  conservation  property 
of  the  reconstruction  (2.3b)  we  get  that 


(3.6a) 


t  e 

=  W-;/‘),E“'*’w>  =  W-;/‘).v>‘>  =  si 
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or  in  matrix  form 


(3.6b) 

It  follows  then  from  (3.5b)  and  (3.6b)  that  for  any  vector  /* 

/*  =  =HRp, 


which  shows  that 


(3.7a) 


HR  =  I, 


and  consequently 


(3.7b)  H(I  -  RH)  =  H-  iHR)H  =  H  -  H  =  0. 


We  turn  now  to  examine  the  error  c*  ^  in  the  prediction  (3.5) 

(3.8a)  =  /*“>  -  /*"*  =  /*"»  -  Rp  =  (/  -  RH)p-^ . 

^From  (3.7b)  it  follows  that 

(3.8b)  /fc*-*  =  0, 

which  shows  that  only  Nk  out  of  the  2Nk  components  of  e*“^  are  independent 
quantities.  In  order  to  get  rid  of  this  redundancy  in  we  introduce  the  Nk  x  2Nk 
matrix  G 

(3-9a)  Gi)  =  (-l)^'’"^a2i-i_>,  (G')a^*x2N* 

which  satisfies 

(3.9b)  HG*  =  0. 
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In  Appendix  A  we  show  that  it  follows  from  the  orthogonality  condition  (2.4c)  that 


(3.10a) 

(3.10b) 

(3.10c) 


HH*  =  GG*  = 
H*H  +  G*G  = 


Using  (3.10b)  and  (3.8b)  we  now  get  that 

(3.11a)  c*-'  =  +  (7*C?)c*-*  =  ^G*(Gc*-')  =  -\G*d\ 

a\*  jap  ap 


where 


(3.11b) 


d*  = 


is  a  vector  of  Nk  components.  Combining  (3.11)  with  (3.8a)  we  get 


(3.12) 


p-1  ^  p-l  ^  ^  ^ 


which  is  the  basis  for  the  following  data  compression  algorithm: 
Given  a  sequence  of  Nq  numbers  u  =  we  set 


(3.13a) 


and  execute 


D.'  for  =  1,2, ...  ,L 


/*  =  Hf^-^ 


(3.13b) 


c‘->  =/*-!- i?/* 
d*  = 


thus  arriving  it  the  multiresolution  representation  of  u 


(3.13c) 
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Starting  from  the  multiresolution  representation  we  recover  u  by  (3.12),  i.e. 


(3.13d) 

(  Do  ioT  k  =  L,L  —  ly. . .  ,  1 
1  /*-'  =  B/*  + 

(3.13e) 

u  =  f. 

The  number  of  quantities  in  the  multiresolution  representation  (3.13c)  is 
No  as  in  the  original  vector  tt  (3.13a).  The  difference  is  that  the  quantities  {df} 
are  expected  to  be  small  in  absolute  value  wherever  the  underlying  function  /(x) 
is  adequately  resolved  on  the  A:-th  grid.  Thus  data  compression  can  be  achieved  by 
setting  to  zero  elements  of  d*'  which  fall  below  some  tolerance  £k-  See  [7]  for  more 
details. 

In  Appendix  A  we  present  the  form  of  the  data  compression  algorithm  (3.13) 
when  we  do  not  assume  the  orthoginality  condition  (2.4c). 

In  the  following  we  present  the  details  for  the  three  cases  that  we  highlight  in 
this  paper. 

Case  1.  Point  values. 

7^(1;/*)  is  the  interpolation  (2.7b).  For  reasons  of  symmetry  we  consider  even 
order  of  accuracy  r  =  2s  and  take  7?.(x;  /*)  in  ,  ar*]  to  be  the  unique  polynomial 
of  degree  (2s  —  1)  that  interpolates  /*  at  the  gridpoints  {x*_, , . . .  ,  x*^,_j }.  In  (3.5) 
we  get  for  1  <  t  <  Nk 


(3.14a) 

/*-■  =  (Rphi  =  ft 

(3.14b) 

ftr-\  = 

t=i 

where 

f  r  =  2  =>01  =  ^ 

^3.14c) 

i 

r  =  4=>/?i  =  $2  = -Jo 
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In  (3.4)  and  (3.9a)  we  get 
(3.15)  Hij  = 

The  multiresolution  representation  (3.13c)  is  obtained  by: 
Set 

(3.16a)  P  =u 


(3.16b) 


'  Do  for  A:  =  1,2, . . .  ,L 

■  4  =  flr-\  -  e;=.  MfK,-,  +  fLiX 


u  is  recovered  from  the  multire, solution  representation  hy 

(  Do  for  A:  =  X,  jL  —  1, . . .  ,1 


(3.16c) 


I  fu'  =  1  <  ■  < 

( + IK,) + <*?  ■  1  <  ■'  < 


(3.16d)  u  =  p 


Case  2.  Cell-averages. 

Using  interpolation  of  order  (2s  +  2)  as  above  for  the  primitive  function  in  (2.11) 
we  obtain  a  reconstruction  of  order  r  =  2s  +  1.  In  (3.5)  we  get  for  1  <  t  <  Nk 


(3.17a)  =  {Rphi-,  =11  +  4 

(3.nb)  /*-■  =  (R/%  =  /*  -  a? 

where 

(3.17c)  zf  ='j2'reifi^t-fi-e) 

/=i 


and 

(3.17d) 


r  =  3  =»  7i  =  -  J 
r  =  5  =►  7l  =  —128  I  ~  128’ 
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note  that  z*  =  0  for  r  =  1. 
In  (3.4)  and  (3.9a)  we  get 


(3.18) 


+  ^2i-lj),  Gij  =  -{62i-i,j  -  62i,j). 


The  multiresolution  representation  (3.15c)  is  obtained  by: 
Set 

(3.19a)  p  =u 


(3.19b) 


’  Do  for  A:  =  1,2, ...  ,£ 

•  ft  = +  ft'%  l<i<N, 

■  4  =  flr-\  -ft-  E:=.  Mfh,  -  ft-,).  1  <  i  <  JVi 


u  is  recovered  from  the  multiresolution  representation  (3.13c)  by 


(3.19c) 


Do  for  A:  =  £,  £  —  1, . . .  ,  1 
Do  for  t  =  1,2,. . .  yNk 


^  Itifi+t  -  fi-t)  +  <^i 

.  fir-\  =  fl+^,  ft' =  ft -A, 


(3.19d) 


u 


Case  S.  Orthogonal  Wavelets 
The  reconstruction  (2.14)  for  the  A:-th  level  is 


(3.20a) 

and  (3.5)  becomes 


Nk 


ft-'  =  W‘)i  =  W-;/*), *>?-')  =  jj^ 
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Using  (3.3c)  and  (2.5b)  we  get  that 


M 

hk-i 


Oti-2j  = 


M 

hk-i 


h:,. 


Recalling  that  hk  =  2hk-i  we  get 


(3.20b) 


R  =  2H\ 


Using  (3.10)  with  (2.5a)  to  express  the  error  in  (3.8a)  we  get  that 


(3.21a)  =  (/  -  RH)p-^  =  (/  -  2H*H)p-^  =  2G*Gp-^  =  2G*d* 

(3.21b)  d^  =  Gp-\ 


The  coefficients  {>/2  o/},  1  <  £  <  2r  of  Daubechies  [4]  are  given  in  the  following  ta¬ 
ble:  Table  1. 


r  =  2 

r  =  3 

r  =  4 

r  =  5 

r  =  6 

'J2a\ 

.482962913145 

.332670552950 

.230377813309 

.160102397974 

.111540743350 

y/2a2 

.836516303738 

.806891509311 

.714856570553 

.603829269797 

.494623890398 

y/2a3 

.224143868042 

.459877502118 

.630880767930 

.724308528438 

.751133908021 

\/2Q!4 

-.129409522551 

-.135011020010 

-.027983769417 

.138428145901 

.315250351709 

v/2a5 

-.085441273882 

-.187034811719 

-.242294887066 

-.226264693965 

y/2a^ 

.035226291882 

.030841381836 

-.032244869585 

-.129766867567 

\/2a7 

.032883011667 

.077571493840 

.097501605587 

y/2a% 

-.010597401785 

-.006241490213 

.027522865530 

\/2a9 

-.012580751999 

-.031582039318 

y/2a\Q 

.003335725285 

.000553842201 

y/2a\i 

.004777257511 

\/2a\2 

-.001077301085 

The  multiresolution  representation  (3.15c)  is  obtained  by; 
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Set 


(3.22a) 


(3.22b) 


/• 


< 


V 


Do  for  k  =  1,2,...  ,L 
Do  for  t  =  1, 2, . . .  ,Nk 


•if = Ej;.(-i)'“ 


u  is  recovered  from  the  multiresolution  representation  by: 


(3.22c) 


(3.22c) 


Do  for  k  =  L,L  —  1,. . .  ,1 
Do  for  t  =  1,2, . . .  ,  JV* 

/2i-l  =  ^Y>t=l[°‘2t-lfi-i  + 

u  =  /°. 


4.  Matrix-vector  multiplication. 

In  this  section  we  describe  a  multiresolution  algorithm  for  the  multiplication  of 
an  iVo -vector  6  by  em  iVo  x  iVo  matrix  A,  which  is  based  on  the  data  compression 
of  j4;  we  denote  the  result  of  this  product  by  the  iVo-vector  c, 

(4.1)  Ab  =  c. 


We  steirt  by  presenting  a  tensor-product  extension  of  the  one- dimensioned  data 
compression  edgorithm  (3.13)  to  the  matrix  case,  in  which  each  column  and  row  of 
the  matrix  are  treated  as  one-dimensional  vectors.  Let  us  set 

(4.2a)  A^  =  A 

and  define  the  Nt  x  Nk  matrix  ^4*  by 

(4.2b)  A^=HA*‘-^H\  jfc  = 
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where  H  is  the  Nk  x  2Nk  matrix  defined  by  (3.4). 

Given  we  form  the  prediction  by 

(4.3a)  A*'-*  =RA^R*, 

where  R  is  the  2Nk  x  Nt  matrix  in  (3.5).  It  follows  from  (4.2b)  and  (3.7a)  that  the 
error  in  this  prediction  E^~^, 

(4.3b)  -  RA^R* 

satisfies 


(4.3c) 


HE^-^H*  =  HA'^-^H*  -  {HR)A^{HRy  =  >1*  -  =  0. 


Consequently,  using  (3.10b)  and  (4.3c)  we  get 


£*-1  =  +  G*G)E'‘-\H*H  +  G*G) 


(4.4a) 


{G*D\G  +  GD^H  +  H*DlG), 


•  rtki 


a 


where  the  Nk  x  Nk  matrices  denote 

(4.4b)  £>*  =  GE’^-^G*,  D\  =  GE^-^H*,  =  HE'^-^G*. 


Thus 


(4.5) 


=  RA^RT  +  ^[G*(Df  G  +  D^H)  +  H*DIG], 


and  we  get  the  following  data  compression  algorithm  for  the  iVo  x  No  matrix  A: 
Set 


(4.6a) 


(4.6b) 


A^  =  A 

Do  for  fc  =  1,2, ...  ,L 

A^  =  HA^-^H* 

< 

Ek-l  =  ^*-1  _ 

,  Df  =  GE*-1G*,  =  GE*‘-'^H\  D*  =  HE^-^G\ 
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The  multiresolution  representation  of  A  is 

(4.7a)  A"*  =  AD])U)]. 

It  is  convenient  to  store  A^^  in  the  form 


(4.7b) 


^2 

I>3 

^2 

dL 

aL 

No 


No 


which  also  shows  that  the  number  of  elements  in  A^^  is  {NoY,  as  in  the  original 
matrix  A  =  A^. 

Starting  from  the  multiresolution  representation  A^^  (4.7),  we  recover  the  orig¬ 
inal  matrix  A  by  (4.5),  i.e. 


(4:8a) 

(4.8b) 


{Do  for  fc  =  L,  L  —  1, . . .  ,1 

=  iL4*i?’  +  |^[G*(D‘G  +  D^H)  +  H*D^G], 

A  =  A^. 


The  elements  of  are  proportional  to  the  local  error  in  predicting  A^~^ 

from  the  fc-th  level  of  resolution  (4.3b).  Therefore  these  elements  are  small  wherever 
the  discretized  function  is  properly  resolved  on  the  jb-th  grid.  Data  compression  can 
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Figure  la 


18 
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be  achieved  by  setting  to  zero  elements  of  which  £ire  smaller  in  absolute 

value  than  some  tolerance  £fc. 

In  Figures  la,b  and  2a, b  we  show  resiilts  of  data  compression  of  two  matrices 
which  Me  the  first  two  examples  in  the  BCR  paper  [l].  In  Figures  la,b  we  show 
the  multiresolution  representation  (4.7b)  of  the  matrix 


(4.9) 


* 

i=j 


with  No  =  512.  The  discretization  in  this  calculation  is  assumed  to  be  by  pointval- 
ues,  i.e.  ff  and  G  are  (3.15)  and  the  reconstruction  is  by  interpolation.  We  take  R 
to  be  (3.14)  with  r  =  6.  Entries  of  which  are  larger  in  absolute  vadue  than 

£k  =  10”^  are  marked  in  black.  The  calculations  in  Figures  la  and  lb  differ  in  the 
treatment  of  boundaries:  In  Figure  la  we  use  periodic  boundary  conditions  while 
in  Figure  lb  we  use  one-sided  interpolation  near  the  boundaries.  The  compression 
rate  (ratio  between  (Nq)^  to  the  number  of  entries  that  are  larger  in  absolute  value 
than  10"^)  is  6.72  for  the  periodic  case  in  Fig.  la  and  8.57  for  the  one-sided  in¬ 
terpolation  at  boundaries  in  Fig.  lb;  the  compression  rate  for  the  wavelet  based 
zdgorithm  in  [1]  is  7.33. 

In  Figures  2a, b  we  repeat  the  ceilculations  of  Figures  la,b  for  the  matrix 
,  lo8l.-No/2[-los|;-Noi21  i  i  /  No/2J  ^  No/2 

(4.10)  Aij  =  \ 

I  0  otherwise. 


Here  the  compression  rates  are  6.11  in  Fig.  2a  and  7.60  for  Fig.  2b;  the  correspond¬ 
ing  BCR  result  is  7.50. 

We  remark  that  the  “BCR  results”  above  are  quoted  from  [1]  in  which  a  different 
normalization  in  (3.3)  is  used.  These  results  show  that  the  BCR  compression  rates 
are  of  the  szune  order  as  the  ones  in  Figures  1  and  2. 

We  turn  now  to  describe  how  to  compute  the  product  Ab  =  c  (4.1)  from  the 
multiresolution  representation  A^^  (4.7b)  of  A.  Multiplying  (4.5)  by  a  vector  6*“^ 
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of  iVfc_i  components  we  get 


+  /f£>5(G6‘-')}, 

from  which  we  see  that  if  for  all  k  we  define 

(4.11b)  6*  = 

(4.11c)  c*  = 

then  (4.11a)  becomes 

(4.12)  =i2c*  +  |^{G*(I>‘(G6*->)  +  I>*(/ffc*-i)]  +  /fI>‘(G6*-i)}. 

It  follows  therefore  that  given  the  (compressed)  multiresolution  representation 
(4.7b)  of  A  we  can  calculate  c  =  Ab  by: 

Set 

(4.13a)  6°  =  b, 

'  Do  for  =  1,2, ...  ,L 
(4.13b)  .  =  = 

.  6‘  = 

evaluate  by  direct  multiplication 

(4.13c)  c^  =  A^b^, 


and  execute 


(4.13d) 


Do  for  A  =  £,  £  —  1, . . .  ,  1 

c*->  =  i2c*  +  ^[G*(Dff*  +  D*s‘)  +  /f(D*<*)], 
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(4.13e) 


c 


c 


0 


Relation  (4.11b)  can  be  thought  of  as  stating  the  proper  scaling  of  the  input 
vector  as  we  go  to  a  coarser  grid.  After  preparing  the  VEilues  of  6*  for  all  the  levels 
(4.13b),  we  start  the  computation  o{  c  =  Ab  by  calculating  its  lowest  resolution 
version  in  (4.13c).  Then  we  proceed  in  (4.13d)  to  successively  upgrade 

c*  by  first  using  the  reconstruction  technique  to  predict  the  value  =  i2c*  for 
the  finer  grid  and  then  correct  this  prediction  wherever  needed  by  the  term  in  the 
curved  brackets  in  the  RHS  of  (4.12). 

If  the  number  of  elements  in  that  are  larger  in  absolute  value  than  the 

tolerance  £k  is  O(iVfc),  and  the  matrices  H,G  and  R  are  banded  (with  constzint 
width),  then  the  number  of  operations  for  each  k  in  (4.13b)  and  (4.13d)  is  0(Nk), 
and  consequently  the  number  of  operations  in  the  multiplication  algorithm  (4.13) 
is  0{Nq). 

It  is  important  to  observe  that  due  to  the  tensor-product  nature  of  this  algorithm, 
the  operations  on  the  rows  axe  independent  of  the  operations  on  the  columns.  This 
enables  us  to  use  Hx ,  Gx  and  Rx  on  the  left  and  different  Hy ,  Gy  and  Ry  on  the 
right.  Modifying  the  relations  (4.2b),  (4.3b)  and  (4.4b)  to  be 


(4.2b)'  A*  =  HxA’^-^h;, 

(4.3b)' 

(4.4b)'  Df  =  GxE^-^G;,  £>2*  =  GxE'‘-'H;,  DI  =  HxE^-^G*y, 


we  now  get  the  following  multiplication  algorithm: 
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Set 


This  extra  freedom  in  algorithm  (4.13)'  can  be  utilized  for  example  to  discretize 
the  integral  transform  (1.1)  by  pointvalues  in  x  and  cell-averages  in  y. 

Next  we  present  details  for  the  three  cases  that  we  highlight  in  this  paper. 

Ca^e  1.  Pointvalues. 

It  follows  form  the  definitions  of  H  and  G  in  (3.15)  that  (4.2b)  and  (4.3b)  become 

(4.14a)  =  l<iJ<Nk, 

(4.14b) 

1  <  ij  <  Nk. 

Using  the  definition  (3.14)  of  R  in  (4.13b)  we  get 

6f  =  =  I*,-'  +  E 1  <  ■  <  N,. 

Algorithm  (4.13)  can  be  expressed  in  this  case  by; 
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Set 


(4.15a)  b°  =  b 

'  Do  for  A:  =  1,2, . . .  ,I« 

(4.15b)  ^  4  =  b^-\  l<i<  Nk, 

■  6?  =  +  EUi  MtUi  +  1  <  *  < 

(4.15c)  =  A^b^, 

Do  for  fc  =  L,  L  —  1, . . .  ,  1 
Do  for  t  =  1, 2, . . .  ,  iVjt 

(4.15d) 

4r-i  =  e;=.  *(-:?+,  -  4-i) + + 02*»*)i 

,  c‘-'  =  c‘  + 

(4.15e)  c  =  c°; 


here  r  =  2^. 

While  writing  this  paper  we  found  out  that  algorithm  (4.15),  edthough  derived 
differently,  had  edready  been  published  in  [2].  Moreover,  it  was  extended  further  in 
[3]  to  integral  transforms  with  zui  oscillatory  kernel  and  to  many-body  problems. 

Case  2.  Cell-averages. 

It  is  convenient  to  introduce  the  operators  fi  and  6, 

(4.16)  fiVi  =  ^(u,  +  Vi_i),  6vi  =  ^(v,  -  Vi_i), 

and  use  the  convention  that,  when  applied  to  two-dimensional  arrays,  superscripts 
X  and  y  denote  operation  on  the  first  and  the  second  index,  respectively.  It  follows 
from  the  definition  of  H  and  G  in  (3.18)  that  (4.2b)  and  (4.3b)  become 


(4.17a) 

(4.17b)  (Df),.>  =  6^6yElj\  (D‘)i,>  =  (D*),.,  = 

1  <  hj  <  Nk. 


26 


Using  the  definition  (3.17)  of  R  in  (4.13b)  we  get 
(4.18)  6*  =  =  2(^6*-'  + 

/=1 

Algorithm  (4.13)  can  be  expressed  in  this  case  by: 

Set 

(4.19a) 

(4.19b) 

(4.19c) 

(4.19b) 

(4.19e) 
here  r  =  2s  +  1. 

Case  S.  Orthogonal  wavelets. 

In  this  case  H  amd  G  are  defined  by  (3.4)  and  (3.9a)  and  the  Daubechies  co¬ 
efficients  (see  Table  1).  Since  R  =  2H*  (3.20b)  Eind  HG*  =  0  (3.9b)  we  get  in 
(4.4b) 

(4.20)  Df  =  =  GA^-^H%  D*  =  HA^-^G*; 


6“  =6 

'  Do  for  fc  =  1,2, ...  ,1/ 

<  sf  =  2/i65-‘ ,  =  28^-^ ,  l<i<Nk, 

^  =  5  •  +  EJ=i  7t{t^+t  -  <*-/)> 

Do  for  A:  =  L,  X  —  1, . . .  ,  1 
Do  for  i  =  1,2, .. .  ,  JVjk 
w  =  c^  +  (D^t^)i 

^  =  E<=i  7<(c?+/  -  cf_,)  -  (I>f 

C2,“\  =  tv  +  2 
>  =  lU  -  2, 
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thus  for  I  <  i,j  <  Nk 


2r  2r 


(4.21a) 

/=1  m=l 

(4.21b) 

/=1  m=l 

(4.21c) 

(^2)i,j  =  E  “”*'^2i-\-/,2>+m’ 

<=1  m=l 

(4.21d) 

(^3)«,i  =  E  (~^)"*Qfni  E  *^^^2»+V2t-l-m- 

m=l  /=! 

Using  R  = 

2H*  in  (4.13b)  we  get  that 

(4.22) 

b‘‘  =  2Hb^-\  s‘‘  =  b^. 

Algorithm  (4,13)  can  be  expressed  in  this  case  by: 
Set 

(4.23a)  6°  =  b, 


/ 


(4.23b)  { 


Do  for  A:  =  1,2, ,  ,1. 
Do  for  t  =  1,2, . . .  ,Nk 


(4.23c)  = 


(4.23d)  < 


(  Do  for  A:  =  L,  L  —  1, . . .  ,  1 
Do  for  t  =  1,2, . , ,  ,Nk 

4-'l  =  2  +  (i>}<‘)i-/l  +  02»(Df  (*  +  X>|»‘)i+a. 

I  4  =  2ELi{“2<(4-c  +  (D3<*)i-«)  -  02/-i(Cf(*  +  i>jH*)i+,), 


(4.23e) 


c  =  c°. 
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Algorithm  (4.23)  is  identical  to  the  BCR  algorithm  (the  “nonstandard  form”)  in 

[I]- 

5.  Stability  and  Efficiency. 

In  this  section  we  examine  the  stability  of  the  data  compression  algorithm  (4.6)- 
(4.8)  and  discuss  the  efficiency  of  the  matrix-vector  multiplication  algorithm  (4.13). 

From  (4.3b)  we  get 

A"  =  E®  -I-  RA^R*  =  -f-  -h  R^A^{R?Y  =  •  •  • 

(5  ll 

=  E°-l-^R’‘E\R'‘y+R^A^(R^y. 

k=i 

Applying  data  compression  to  A^^  (4.7)  we  get  truncated  matrices  E^  which  result 
in  A®  in  (5.1).  Denoting 

(5.2a) 

we  thus  get 

t-i 

(5.2b)  A° -A°  +  '^R^S'‘iR’^y, 

fc=i 

which  shows  that  each  column  and  row  in  £*  are  amplified  by  R*.  For  discretization 
in  [0, 1]  R*  in  (5.1)-(5.2)  should  be  interpreted  as 

(5.3a)  R*  =  R,  •  R2  •  •  •  Rib 

where  Rm  is  the  2Nfn  x  Nm  matrix  in  (3.5);  for  discretization  in  (— oo,  oo)  R  is  an 
infinite  matrix  and  R*  should  be  interpreted  as  the  k-ih.  power  of  R,  i.e. 

(5.3b)  R*  =  (R)*. 

Let  c  denote  the  imit  sequence  corresponding  to  a  partition  of  the  real  line  into 
intervals  of  size  1  with  integer  endpoints, 

c/  =  —<x>  <  £  <  oo, 
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and  consider  successive  applications  i2*e,  k  — »  oo.  For  exaimple  when  R  is  the 
piecewise  linear  interpolation  (3.14)  with  r  =  2  we  get 


Table  2. 

X  =  —1  X  =  0  X  =  1 

c  0  1  0 

iZc  0  I  1  i  0  0 

iZ^cOOjfflfJJOOO 

0  0  0  i  I  I  I  I  I  I  1  I  I  f  I  I  t  J  0  0  0  0  0 

Clearly  here 

(5.4a)  =  ,(2'*,), 

f  1  -  kl  kl  <  1 

(5.4b)  ti(x)  =  < 

(  0  otherwise 

We  observe  that  r]{x),  the  “hat  function”,  is  the  solution  of  the  dilation  equation 
(5.4c)  Ti{x)  =\ri{2x  -  1)  +  ri{2x)-\-\r){2x  +  1), 

the  coefficients  of  which  are  given  by  at  =  (i2e)<. 

The  limiting  process  i2*e,  k  —*  oo,  has  been  studied  by  Deslauriers  and  Duboc 
[5]  and  Dyn,  Gregory  and  Levin  [6]  for  interpolating  R,  and  by  Daubechies  [4]  for 
orthonormal  wavelets,  R  =  2H*  (3.20b).  As  in  the  example  above  they  found  that 
the  hmiting  process  is  convergent  in  the  sense  that 

OO 

(5.5a)  lim  ^  (i2*c)jA'jo,i)(2*x  -  j)  =  »7(x), 

j=-oo 

where  A'[o,i)  is  the  characteristic  fimction  of  [0, 1)  and  the  convergence  is  uniform 
in  X.  The  limit  f7(x)  is  a  continuous  function  of  compact  support  which  satisfies  the 
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dilation  equation 


(5.5b) 

'li^)  -  Yi  =  (^e)/, 

j 

and 

€ 

(5.5c) 

^-‘  =  (,,2*,.(2‘--3))=(fiMr 

Since  t7(x)  is  continuous  and  of  compact  support  it  follows  from  (5.5a)  that 

OO 

(5.6a)  supjfc{2~*  ^  l(i?*e)j|}  <  const ,  supy,*|(i?*e)j|  <  const. 

i=-oo 

Consequently  we  get  for  the  matrix  norms 


(5.6b) 


||i2‘||oc  <  Coo,  ||i2*lli  <2*-Ci. 


We  return  now  to  the  stability  analysis  (5.2)  of  the  data  compression  algorithm. 
Setting  to  zero  elements  of  ZJ*  (4.7)  which  fail  below  the  tolerance  e*  we  get 

(5.7a)  <  const  •  £*+1, 

(5.7b)  11^*11;,  <C-iV*.£fc+i,  p=l,oo. 

For  each  term  in  (5.2b)  we  now  get  for  both  the  Li  and  Loo  norms  that 

(5.8a)  ||JJ‘£‘(R*)-||  <  ||il‘|U||iJ*||i  -C-Nf  £j+,  <  CC^C,  •  JV„  •  £t+, 

=  C  •  No  •  €k+i, 

and  consequently 

(5.8b)  J;££; 

"»  h.i 

this  shows  the  stabihty  of  the  data  compression  algorithm. 

In  the  numerical  experiments  shown  in  Figures  la,b  and  2a, b  for  the  matrices 

(4.9)  and  (4.10)  we  have  used  Ck  =  e  =  10"^  (here  ho  =  1)  and  computed 

(5.9)  i;,(£)  =  IKA”  -  A)5||,/||6||,.  p=l,oo 


31 


for  a  randomly  generated  vector  6;  for  purposes  of  compsirison  we  also  computed 
C'p{0)  which  corresponds  to  running  the  program  with  £  =  0  and  thus  shows  the 
effect  of  round-off  error.  In  Table  3  we  show  the  results  for  the  case  where  R  is  the 
interpolation  (3.14)  with  r  =  6. 


Table  3. 


case 

boundary 

ratio 

^i(io-") 

f'oo(10  '^) 

i>i(0) 

i^oc(O) 

(4.9) 

periodic 

6.72 

one-sided 

8.57 

7.52  X  10-® 

4.41  X  10-5 

2.77  X  10-5 

(4.10) 

periodic 

6.11 

1.62  X  10-« 

1.82  X  10-5 

4.76  X  10-® 

9.15  X  10-® 

one-sided 

7.60 

1.46  X  10"® 

It  seems  to  us  that  the  convergence  of  R^e  to  a  continuous  r]{x)  stems  from  the 
conservation  property  HR  =  I  and  the  accuracy  requirement  (2.3a)  with  r  >  2. 
Therefore  we  expect  the  reconstruction  from  cell-averages  (3.17)  to  also  satisfy  the 
relations  (5.5),  (5.6).  In  Appendix  B  we  prove  convergence  of  the  limiting  process 
(5.5a)  for  reconstruction  from  cell-averages  under  the  assumption  that  the  corre¬ 
sponding  limiting  function  for  the  interpolation  in  (2.11)  is  continuously  differen¬ 
tiable. 

In  Table  4  we  repeat  the  calculation  of  Table  3  for  the  reconstruction  from  cell- 
averages  (3.17)  with  r  =  5. 
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Table  4. 


case 

boundzury 

ratio 

i>i(10-'') 

Mio-'') 

i>i(0) 

i>oc(0) 

(4.9) 

periodic 

5.71 

6.03  X  lO""^ 

7.83  X  lO-"^ 

5.66  X  lO-"^ 

4.39  X  lO-"^ 

one-sided 

6.71 

1.03  X  10"® 

1.55  X  10-« 

9.87  X  lO""^ 

9.52  X  10-’ 

(4.10) 

periodic 

6.29 

4.00  X  10-^ 

5.97  X  10"^ 

3.50  X  10-^ 

3.06  X  10-^ 

one-sided 

7.53 

2.76  X  lO-"^ 

6.09  X  10"'^ 

1.73  X  lO-"^ 

3.06  X  lO-"^ 

We  turn  now  to  discuss  the  question  of  efficiency.  If  a(i,  y)  is  a  function  that  has 
isolated  regions  of  large  variation  then  its  discretization  on  a  uniform  grid  results  in 
a  matrix  A  which  is  actually  over-resolved  in  most  of  the  computational  domain.  In 
this  case  it  pays  to  use  multiresolution  algorithms  as  they  offer  the  efficiency  of  an 
adaptive  grid  method  without  the  complicated  logics  that  is  associated  with  such  a 
calculation.  In  applying  multiresolution  algorithms  to  matrix- vector  multiplication 
there  is  another  important  consideration:  The  computationed  effort  of  prepairing  the 
representation  (4.7)  may  be  greater  than  a  direct  application  of  the  matrix 

i4  to  a  single  input  vector  b.  Therefore  it  makes  sense  to  use  algorithm  (4.13)  only 
when  the  computationzil  task  calls  for  an  application  of  the  same  matrix  to  many 
input  vectors  and/or  there  is  apriori  knowledge  of  the  location  of  regions  of  large 
variation. 

An  important  class  of  applications  is  the  czdculation  of  integred  transforms  (1.1) 

(5.10)  u(x)=  f  K{x,y)v{y)dy, 

Jo 

where  the  kernel  K{x,  y)  is  smooth  except  for  curves  y»{x)  at  which  it  has  integrable 
singularity.  To  each  grid  of  size  hk  we  associate  a  finite-dimensional  approximation 
K^{x,y)  to  the  kernel  K(x,y) 


(5.11a) 

(5.11b) 


Nk  N„ 
i=l  j=l 

Jo  Jo 


33 


and  a  finite-dimensional  approximation  «*(x)  to  the  output  u(i) 

N*  Nk 

(5.12a)  «*(i)=  /  K*‘{x,y)v{y)dy  = 

»=i  >=i 

(5.12b)  v(y)r}^{y)dy. 

Here  r}(x)  is  the  limiting  function  in  (5.5)  and 

(513a)  =  f}(— -  ey 

^From  (5.5c)  with  A:  =  0  we  get  that 

(5.13b)  jT}{x)<p{x-j)dx  =  6oj 

and  thus  by  scaling 

(513c) 

Using  (5.13c)  in  (5.12a)  we  get 

Nk 

(5.14a)  uf  =  /H  ^  l<i<  Nk, 

j=i 

(5.14b)  0?  =  (u‘,<^f). 

Applying  the  data  compression  algorithm  to  the  matrix  and  setting  to  zero 
elements  of  {I?*}  that  fall  below  et  we  get  from  (5.8b)  that 

L 

(5.15a)  <C5];efc. 

It  follows  therefore  that 

(5.15b)  ||a"  -  u»||  =  ||/.„(K«  -  K-")v“ll  <  C(52ej)||S»||. 

*=1 


34 


The  prediction  error  E^j  ^  (4.3b)  cein  be  estimated  by 

(5.16a)  |£5-  I  < 

where 


(5.16b) 


_ 


=  (lS  +  l?fl).i- 


dx^ 


dy^ 


If  the  kernel  K{x^y)  is  such  that  at  a  distamce  A  from  the  singularity  can  be 
bounded  by 


(5.17a) 


Ck 

A*" 


then  it  follows  from  (5.16)  that  except  for  a  band  of  width  B  aa*ound  the  singularity, 
all  the  elements  of  satisfy 

(5.17b)  |Bj-|  < 

Choosing  =  £  in  (5.15)  amd 

(5.18a)  B  = 

we  get  that  the  number  of  nonzero  elements  in  the  compressed  (4.7)  is 

proportioned  to  2BNk,  and  the  cumulative  error  (5.15b)  is 


(5.18b  l|u°  -  u°||  <  C||t;®|l£log2  Nq. 

This  error  can  be  made  arbitrarily  small  by  taking  a  wider  band  B  in  (5.18a). 

Beylkin,  Coifman  and  Rokhlin  [1]  point  out  that  the  estimate  (5.17a)  is  satisfied 
by  the  kernels  of  Calderon- Zygmund  operators  and  pseudo-differential  operators. 
In  this  case,  taking  into  account  the  actual  decay  of  the  prediction  error  away  from 
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the  singulairity  to  sharpen  the  estimates  in  (5.7),  the  log2  Nq  factor  in  the  RHS  of 
(5.18b)  C2in  be  removed. 

6.  Summary  and  Conclusions. 

In  this  paper  we  have  presented  a  class  of  multiresolution  algorithms  for  data 
compression  and  matrix-vector  multiplication.  In  constructing  this  class  we  have 
introduced  subclasses  of  different  discretizations.  Each  subclass  corresponds  to  a 
particuleu’  choice  of  <p(x)  in  (2.2);  tp(x)  is  assumed  to  be  a  solution  of  a  dilation 
equation  and  to  satisfy  the  orthogonality  condition  (2.4c).  Members  of  each  sub¬ 
class  of  discretization  correspond  to  different  reconstruction  procedures  7l(x;  /);  the 
reconstruction  is  assumed  to  be  conservative  (2.3b)  and  to  depend  linearly  on  the 
discrete  data  /. 

We  have  paid  special  attention  to  the  subclasses  of  discretization  corresponding 
to  pointvalues  and  cell-averages  because  of  their  simplicity.  The  wavelet  based  al¬ 
gorithms  [1]  are  also  included  in  this  class  but  in  a  “diagonal”  fashion:  In  each 
subclass  of  discretization  corresponding  to  a  ip(x)  which  satisfies  the  moment  con¬ 
dition  (2.15),  there  is  a  wavelet  based  algorithm  corresponding  to  the  reconstruction 
iZ  s=  2ff*  (3.20b).  For  example  the  wavelet  based  algorithm  for  r  =  1  (Haar  basis) 
is  in  the  subclass  of  cell- averages. 

The  rate  of  compression  and  the  stability  properties  are  about  the  same  for  all 
algorithms  of  this  class  with  the  seime  order  of  accuracy.  What  matters  therefore  in 
choosing  an  eJgorithm  is  simplicity,  operational  count  and  suitability  to  the  pau-ticu- 
lair  application;  under  simplicity  we  also  include  heindling  of  boundaries.  Comparing 
wavelet  bzised  algorithms  to  those  of  pointvalues  and  cell-averages  of  the  S2une  or¬ 
der  of  accuracy  r,  we  find  the  wavelet  based  algorithm  to  be  considerably  more 
expensive  because  of  the  larger  support  (2r)  and  lack  of  symmetry  and  that  the 
handling  of  boundaries  is  not  as  simple.  In  comparing  cell-averages  to  pointvalues 
we  find  cell-averages  to  be  more  suitable  for  discretization  of  kernels  with  integrable 
singularity. 
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Appendix  A. 


Let  P  denote  the  symmetric  matrix 


(Ala) 


P  =  H*H  +  G*G 


where  H  cind  G  are  (3.4)  and  (3.9)  respectively.  A  direct  calculation  shows  that 


(A.16) 


Pij  =P(I*  -  Jl) 


where  for  m  integer 


(A.lc) 


I 


p(2m  —  1)  =  0 

p(2m)  =  ^ojkafc+2m . 
k 


Let  us  2issume  now  that  P  is  an  invertible  matrix.  It  follows  then  from  (3.8b) 
that 


(A.2a)  e*-'  =  P'^Pe^'^  =  +  G*G)e*'‘  =  p-^G’Gc*"* 

=  P'^G’d* 


where 


(A.26) 


d  ‘  =  Gt 
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Replacing  relation  (3.11)  by  (A. 2)  we  get  that  the  encoding  part  (3.13b)  of  the 
data  compression  algorithm  (3.13)  remains  the  same,  but  the  decoding  part  (3.13d) 
becomes 


(A.3) 


Do  for  k  =  L,L  —  1, ..,  1 
7*"^  =  ii7*  +  p-^c?*d*. 


The  orthogonality  condition  (2.4c)  implies  that 


(A.4) 


P  =  p(0)/  =  lap/ 


which  brings  us  back  to  (3.13d). 

As  an  example  for  the  nonorthogonal  csise  let  us  consider  the  “hat  function”  (p(x) 


r  l-|x|  0<  la:|  <  1 
(A.5a)  (p(x)  =  < 

(  0  otherwise 

v/hich  satisfies  the  dilation  equation 


(A.5b)  (p(x)  =  ^[(p{2x  -  1)  +  2(p(2x)  +  <p{2x  +  1)]. 

In  this  case  the  only  nonzero  elements  of  P  eire 


(A. 6)  Pi,»  —  g  ,  Pi,i±2  —  • 

Thus  P  is  diagonally  dominant  and  hence  invertible. 
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Appendix  B. 


In  this  appendix  we  use  the  interpolation  results  of  [5]  and  [6]  in  order  to  prove 
convergence  of  the  limiting  process  (5.5a)  for  cell-averages  with  the  symmetric  re¬ 
construction  (3.17). 

Let  R  denote  the  matrix  (3.5)  corresponding  to  the  central  interpolation  (3.14) 
and  let  fi{x)  denote  the  limit  function  in  (5.5a).  ^(x)  has  its  support  in  |x|  <  r  —  1 
where  r  is  the  order  of  accuracy  of  the  interpolation.  For  r  =  2,  fi{x)  is  the  “hat- 
function”  (5.4b)  which  is  only  Lipschitz-continuous;  for  r  =  4, 6,  ^(x)  is  continuously 
differentiable. 

Let  5”*  denote  the  “step-sequence” 

(  0  j  <  m  —  1 

57*  = 

I  1  J  >  ’Tl. 

The  limiting  process  corresponding  to  is  also  convergent  and  we  denote  its 

limit  by  C(^)-  Since 

e  =  5®  -  5^ 

we  get  that 

(B.la)  fj{x)  =  C(x)  -  C(a:  -  1) . 

It  is  easy  to  see  that 

'0  X  <  — r  -f- 1 

(B.lb)  C(^)  =  -  E/=o^^(^-0  -r-l-l<x<r-2 

.1  r  -  2  <  X 

and  thus  C(x)  has  at  least  the  same  smoothness  as  fj{x). 
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We  turn  now  to  express  the  limiting  process  R^e  for  the  reconstruction  from 
cell-averages  (3.17)  in  terms  of  CC^c).  FVom  (5.5c)  and  (2.11)  we  get  that 


(B.2a) 


f^jt  .  C0'2-‘)  -  C((i  - 1)2-*) 

(«  '))  = - ^ - • 


Since  C'{x)  is  continuous  and  of  compzict  support  we  get  that 


(B.26)  r)(x)  =  ^1^  ^(B‘e)j  Xlo-i)2-*,  i2-*](a:)  =  C'(x) 

j 

and  that  the  convergence  is  uniform  in  x.  From  (B.lb)  it  follows  that  t/(x)  has  its 
support  in  — r  -I- 1  <  X  <  r  —  2;  from  (B.la)  and  (B.2b)  we  get  that  ri{x)  is  related 
to  ^(j)  by 


(B.3)  ^'(x)  =  t;(i)-t7(x-1). 

We  remark  that  for  r  =  2  we  get  for  all  k  that 

X((i-i)2-*,  >2-*](a;)  = 
j 

where  (p{x)  is  the  “box-function”  (2.8a)  (note  that  the  order  of  accuracy  of  the 
reconstruction  from  the  cell  averages  is  r  —  l).  Thus  rj{x)  —  v’(x)  and  we  get  formal 
pointwise  convergence  of  (B.2b)  although  ri{x)  is  discontinuous. 
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